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Abstract. - Inspired on the continued-fraction technique to solve the classical Fokker-Planck 
equation, we develop continued-fraction methods to solve quantum master equations in phase 
space (Wigner representation of the density matrix). The approach allows to study several 
classes of nonlinear quantum systems subjected to environmental effects (fluctuations and dis- 
sipation), with the only limitations that the starting master equations may have. We illustrate 
the method with the canonical problem of quantum Brownian motion in periodic potentials. 



Introduction. The phase-space formulation of quantum mechanics [1] has received a 

renewed attention in the last decades because it allows to employ notions and tools of classical 
physics in the quantum realm. The central object in this formulation is the Wigner function 
W(x,p), which is a phase-space representation of the density matrix p(x,x'). In a closed 
system, the equation governing its time evolution is equivalent to the Schrodinger equation, 
while in the classical limit it reduces to the Liouville equation for the phase-space distribution. 

In an open system, i.e., in contact with the surrounding medium, the system experiences 
generically dissipation, fluctuations, and decoherence [2]. Thus, the study of quantum open 
systems is of interest in many areas of physics and chemistry. Typically one is interested in 
a system consisting of a few relevant degrees of freedom coupled to its environment, or bath, 
which has a very large number of degrees of freedom. 

Under certain conditions, the dynamics of an open system can be formulated in terms of a 
quantum master equation for the (reduced) density matrix [3,4]. In the Wigner representation 
the master equation goes in the classical limit over the Klein-Kramers equation (Fokker- 
Planck equation in phase space [5] ) . Recall that phase-space dynamics also includes problems 
reducible to "mechanical" analogues: Josephson junctions, certain electrical circuits, chemical 
reactions, etc. The most powerful technique to solve the Klein-Kramers equation is the 
continued-fraction method. The distribution is expanded in a basis and the equations for the 
expansion coefficients derived (moment equations). These equations, the so-called Brinkman 
hierarchy, have a tridiagonal structure that can be exploited to find the solution with continued 
fractions [5] . The method has also been successfully extended to rotational Brownian motion 
problems of classical spins and dipoles (see, for instance, [6] and references therein). 
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To deal with quantum dissipative systems, however, is a more complex task. Quantum 
master equations, except in some simple cases, are difficult to solve. Exact expressions for 
the time evolution of the density matrix obtained by path- integral methods are available [7] . 
However, they are in general difficult to evaluate, even numerically, because the propagating 
function is highly oscillatory. Finally, quantum Monte Carlo simulations can in principle be 
used, but they are computationally complex and suffer from the (dynamical) sign problem. 

This situation constitutes a strong motivation for the development of alternative methods 
for quantum dissipative systems, even if they are valid only for specific classes of problems. 
Restricting our attention to the quantum master equation approach, Shibata and co-workers 
[8] developed continued-fraction methods for a quantum spin in contact with a thermal bath. 
Vogel and Risken [9] employed similar techniques in quantum nonlinear optics. However, to 
the best of our knowledge, this method has not been tried for the quantum master equations 
in genuine phase-space problems. 

In this article we derive the hierarchy of moment equations corresponding to a generic 
master equation for the Wigner function of the Caldeira-Leggett type [10]. Then we show 
that, in various cases of interest, the continued- fraction method for the classical problem 
can be adapted to solve this quantum hierarchy, without a sizable increase in sophistication 
or computational complexity. This results in a suitable nonperturbative technique to study 
nonlinear quantum systems subjected to dissipation and thermal fluctuations. Finally, we 
illustrate the method with the problem of quantum Brownian motion in periodic potentials. 

Quantum master equations in phase space. We shall consider the following generic 

master equation for the Wigner function of a particle of mass M in a potential V(x) [10, 11]: 



d t W = 
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The first two terms (Poisson bracket) generate the classical reversible evolution (Liouville 
equation). The next accounts for irreversible effects due to the coupling to the environment 
(7 and D pp are the dissipation and diffusion coefficients). The mixed diffusion term D xp d 2 W 
is heuristically related with the colour of the quantum noise. Finally, the infinite series of 
derivatives dp W gives the quantum contribution to the unitary evolution of the closed 
system. Conditions under which this series can be truncated are sometimes discussed, implic- 
itly assuming that solving an infinite-order partial differential equation would be a hopeless 
task. We shall however keep these terms, which can be specially important in nonlinear sys- 
tems [12]. To conclude, in the phase-space formulation, the quantum- mechanical average of 
any operator A is simply obtained from the corresponding classical variable A(x,p) as 

(A) = Tr(pi) = jdxdpW{x,p)A(x,p) , (2) 

i.e., multiplying A(x,p) by the "distribution" W(x,p) and integrating over the phase space. 
Conditions of validity for the description in terms of quantum master equations are dis- 
cussed in detail in [3,4]. We only mention that equations of the type of eq. Q} are usually 
derived under assumptions like semiclassical (high-temperature) bath, weak system-bath cou- 
pling, etc. Nevertheless, the quantum master equation recently derived for strong coupling, 
and valid for all hj/ksT [13], involves terms structurally similar to eq. Q (with x-dependent 
coefficients). In addition, the phase-space representation of the celebrated Lindblad master 
equation [14] is obtained by simply adding to eq. 0J terms of the form d x (xW) and d 2 W. 
All these types of terms can be readily incorporated in the treatment below. 
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We shall use scaled variables based on a reference length Xo and a frequency Hq (e.g., 
energy scaled by E — MQqXq, action by Sq = MQ Xq, etc.). Besides, D pp (= ^Mk^T in 
a high-T bath) is handled as an effective temperature k^T c a — D pp /jM and p is rescaled 
by (Mfc B Toff) 1/2 . Then V(x) enters divided by_fc B T off and 7 as j T = 7(Ma;§/fcBT cff ) 1 / 2 . 
Finally h is introduced as h/So = 2tt/K, being K related to the ordinary Kondo parameter 
by K = (-f/£lo)K. AdB = tt'Jt/K is the thermal de Broglie wave length (in units of Xo)- 

Derivation of the quantum hierarchy of moment equations. - In the expansion into com- 
plete sets approach to solve kinetic equations in non-equilibrium statistical mechanics [15, 
p. 175], the distribution is expanded in an orthonormal basis {ipn} and the equations of mo- 
tion for the coefficients C n derived. In our case we expand W into Hermite functions ip n (p) 

w { x, P) = w ^c n{ x)M P ), -=g^e-* ( ^ *. = e ^00jr ■ < 3 > 

The Boltzmann-like factor w is extracted for mathematical convenience and involves the auxil- 
iary parameter < 77 < 1/2 and the effective potential Q(x) [usually a scaled version of T^(x)]. 
The Hermite basis has a number of advantages [5] ; not the least is the handling of the deriva- 
tives dp' by means of the associated "creation" and "annihilation" operators b + — — d p + ip 
and b — dp + ^p, which obey the ordinary boson commutation rule \b , b + ] = 1 

From the orthonormality of the {ip n } we have C n (x,t) — Jdpip n (W/w). Differentiating 
with respect to t and using the master equation d t W = C W, one gets the dynamical equations 
for the "coefficients" in the form d t C n = J2 m QnmC m . To get the (operator) matrix elements 
Qnm = ^ dp xp n {w~ l C w)xp m one takes advantage of results for normal ordering of (b + ±b) s 
[16]. This leads to the following quantum (Brinkman) hierarchy 

[(n-l)/2] 
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+ y/{n +l)(n + 2) 7+ C n+2 + ]T K+ 2s+1) V^ +1 )] C n+(2s+1) . (4) 
The auxiliary parameters introduced are combinations of 7t, t], f]± = 77T 1/2, and a = ?/-?/+: 

7± = It »?±(1 - V±) , 7o = It [(??+ -0-) + 2n(?y - a)] . (5) 

The operators on the x dependence read [the T's act only on V(x) not on the C m (jc)'s] 

D ± = d± (8 X - *') , d± = 1 + r]±D xp , A = o-\ 2 dB 3 2 x , (6) 



r^ = 4^«We-^V 1 (-m,2 S + 2;A) J «W = ^^ ^ ^ , (7) 

being i-Fi(a, c ; z) the confluent hypergeometric (Kummer) function and m = n — (2s + 1). 

In the classical case only the s — terms survive in the above sums; if we set the usual 
77 = 1/2, eq. (0J reduces to the 3-term Brinkman hierarchy associated to the classical Klein- 
Kramers equation [5]. However, except for polynomial potentials V^ = 0, Vs > S (harmonic 
oscillator, Duffing oscillator, etc.), the quantum terms lead to an infinite range of coupling in 
the index n. This prevents the use of the custom continued- fraction method (based on the 
n-recurrence) . In the following we shall exploit the expansion in the position basis to show 
how this problem can be circumvented. 
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Expansion into the position basis. The coefficients C n are still functions of x. By 

expanding them in an orthonormal basis {u a (x)} (to be specified for a given problem), we get 
a two index recurrence involving matrix elements of the form A a p = Jdx u* a A(x, d x ) up 



C n (x) = Y,C%u a (x) - -CI = J2H(Q™)«b C ™ 



(8) 



m /3 



Note that the apparent discreteness of this equation is not a consequence of any approximation, 
e.g., grid discretisation. 

Now, introducing appropriate vectors and matrices, the 2-index scalar recurrence can be 
cast into a vector recurrence in the x-basis index (TV is a large truncation index of the p basis) 
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To appreciate better the structure of the "coefficients" Q a p, their general expression is de- 
composed into the "free" and potential contributions Q Q/ 3 = Q f a a + Q^a- For Q { a a we find 
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while the part due to V(x) (including the quantum terms) has the alternate dense structure: 
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The important point is that, for a given problem, with an appropriate choice of the basis 
functions {u a (x)}, the matrix elements in the indexes (a, (3) may vanish when the second 
index lays at a certain "distance" of the first. Then, the recurrence relation J3J has a finite 
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Fig. 1 - Wigner functions in the absence of dissipation for a particle in a sinusoidal potential (one 
period and the two adjacent half-periods are displayed). The reduced Kondo parameter K = 2-K^So/h) 
is 1, 1.57, and 2.22 (left to right), corresponding to q = 1/10, 1/4, and 1/2 of [17, fig. 3]; q = (K/ty) 2 . 

coupling range, e.g., c a — Q Q , Q -iC Q _i +Q a , a c a +Q a ,a+ic a +i, and can always be tackled by 
(matrix) continued- fraction methods [5]. Once the equations are solved, we can reconstruct 
the Wigner function W(x,p) from the expansion coefficients C" and get any observable from 
the general expression J2J for the averages. 

Application to quantum Brownian motion in periodic potentials. - In this problem one 
can anticipate that the coupling in a will coincide with the number of harmonics in V{x): 1 
for cosine potential (3-term recurrence), 2 for a simple ratchet potential (5-term recurrence), 
etc. Let us introduce the Fourier expansion of V'(x) and plane waves as basis functions: 



V\x) 



V'e il 
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(10) 



It is convenient to extract the external force F from V(x) [including it in D ± = d±(d x — $') — 
n±F] and set the auxiliary potential <& proportional to the periodic part <& = eV . Then, to 
obtain all matrix elements in Q Q/ 3 we only need (d x ) a p = iaSafi an d [V"^*^] aj g = V^Jg, whence 
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Here q = a-(3 and G s n {z) = \n ( n ] {q)\ exp(-z/2) iFi(-m, 2s+2 ; z) [cf. eq. Q], being n { n\q) the 
previous quantum coefficient with AdB replaced by qXdB- Note that a = ??-?/+ = rj 2 — 1/4 < 0. 
Then the argument of G s n is positive and exp(— z/2) acts as a Gaussian regularisation factor 
[this is one advantage of the inclusion of the parameter rj in eq. IJ3J)]. This factor reduces 
the weight of the off-diagonal terms inside Q Q| a, contributing to the numerical stability and 
allowing to study the deep quantum regime (AdB ^ ^o)- 

Let us give some examples. Kandemir [17] obtained the eigenfunctions of the Schrodinger 
problem of a quantum particle in a sinusoidal potential in terms of Mathieu functions. They 
were represented with the Wigner function associated to the density matrix of the pure states 
p(x,x') — \&(x)\&*(x'). Setting the damping close to zero (7/fio = 10 -6 ) we reobtain his 
results for the ground state (fig.QJ. It is nicely seen the evolution from delocalised solutions for 
the lowest K to localised ones as K oc So/H increases (the system becomes more "classical"). 

After checking the connection with the Hamiltonian limit, let us include dissipation and 
temperature. Chen and Lebowitz [18] studied the transport properties of underdamped parti- 
cles in a cosine potential. For low forces they obtained a free-particle like behaviour (p) oc F/j. 
Increasing F, the wave-vector associated to p reaches the first zone boundary, where Bragg 
scattering reduces the velocity. Eventually, for larger forces, p's corresponding to states in the 
next band become available, and the free-particle behaviour is progressively recovered. 
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Fig. 2 - Left panel: (p) vs. force for underdamped particles at T = 0.5. The reduced Kondo parameter 
is K = 0.7854, which corresponds to Q q = 4 of [18]; fl q = ■n/K. Right panel: Side view of the Wigner 
function (W vs. p, Vx) for 7 = 0.01 at F/7 = 4.5, 4.75, 5, 5.25, and 5.5 (left to right). 



Since high T/7 approximations are involved in [18], we set D pp = jMk^T and D xp = in 
the master equation Q. Solving it with the continued- fraction method we reobtain the effect 
described above (see fig- EJ - However, with this method we also get the Wigner distribution. 
Representing the obtained W vs. p in the range where the curves start to rise again (.F/7 ~ 5), 
we find two peaks at values associated roughly to the lowest bands (p ~ 2 and 6), supporting 
their interpretation. In analogy with the classical case, we would speak of multistability 
between the "locked" (p ~ 0) and two quantum "running" solutions. In the classical case the 
running solutions have a wiggling structure along x [5, fig. 11.22] (the particle slows down near 
the potential maxima) . Displaying the side view of W , instead of the true marginal distribution 
P(p) — Jdx W(x,p), we can appreciate the straight structure (substrate insensitive) of the 
quantum running solutions. 

Before concluding let us consider an example of dynamical response. Figure shows 
the linear dynamical susceptibility vs. frequency of the ac field of a particle in a sinusoidal 
potential. In the classical case, the line-shape x"(w) broadens and extends to w's lower than 
the frequency of oscillation near the bottom of the potential wells (ui = 1 in scaled units). The 
reason is the dependence of the oscillation period on the amplitude of oscillation in anharmonic 
potentials, which yields a non-dissipative contribution to x"( w ) [5]- 

In the quantum regime (approached by decreasing the Kondo parameter K oc jK) the 
x"(w) curves develop a multipeaked structure, which could be loosely interpreted as the 
nonlinear oscillations becoming quantised. We have calculated analytically the effect of the 
anharmonic terms of a cosine potential (those beyond ex a; 2 ) on the energy levels of the har- 
monic oscillator part. Doing this with first-order perturbation theory, we find the dependence 
on the energy-level index k acquired by the distance between levels AEk+i.k- The correspond- 
ing frequencies AEk+i,k/h agree quite well with the location of the main peaks (fig. El right 
panel), substantiating the above interpretation. Finally, when the temperature is lowered the 
thermal population of the higher levels is reduced. Then, the transitions between these levels 
("large amplitude oscillations") become less probable and the corresponding peaks reduced. 

Discussion. - We have adapted the continued-fraction method of classical Brownian mo- 
tion to tackle quantum master equations in phase space. As the classical counterpart, the 
method is limited to a few degrees of freedom. Besides, it is more problem dependent, as it 
uses the x-recurrence and requires a good position basis. Nevertheless, it inherits most advan- 
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Fig. 3 - Imaginary part of the dynamical susceptibility. Left panel: T = 0.05 and K = 200 with 
7 = 0.01, , 0.003, 0.001, 0.0003, 0.0001 (bottom to top). Right panel: 7 = 0.0001 with T = 0.05 and 
T = 0.025 (lower curve). Vertical lines: transition frequencies of the associated anharmonic oscillator. 



tages from the classical method: efficiency, accuracy, and the obtaining of the "distribution" 
W(x,p), from which any observable can be computed. Besides, it is essentially nonperturba- 
tive and specially apt for stationary solutions (static and dynamic). The eigenvalue spectrum 
is not required, which is desirable when dealing with non-bounded Hamiltonians, continuous 
spectra, etc. The connection with the classical results is always possible and attained in a natu- 
ral way. Indeed, the knowledge of the classical phase-space structure (separatrices, attractors, 
etc.) plus the visualization of the Wigner function, can provide valuable insight into complex 
problems involving nonlinear quantum systems subjected to dissipation and fluctuations. 

* * * 

This work was supported by DGES (Spain), project BFM2002-00113. Valuable discussions 
with F. Falo and D. Zueco are gratefully acknowledged. 

REFERENCES 



[1] Hillery M. et al, Phys. Rep., 106 (1984) 121. 

[2] Weiss U., Quantum Dissipative Systems (World Scientific, Singapore) 1993. 

[3] Karrlein R. and Grabert H., Phys. Rev. E, 55 (1997) 153. 

[4] Ankerhold J., in Irreversible quantum dynamics, edited by F. Benatti and R. Floreanini 
(Springer, New York) 2003, Vol. 622 Lecture notes in physics, p. 165, cond-mat/0309284 

[5] Risken H., The Fokker-Planck Equation (Springer, Berlin) 1989. 

[6] Kalmykov Y. P. and Coffey W. T., Phys. Rev. B, 56 (1997) 3325. 

[7] Grabert H., Schramm P., and Ingold G.-L., Phys. Rep., 168 (1988) 115. 

[8] Shibata F., J. Phys. Soc. Japan, 49 (1980) 15. 

[9] Vogel K. and Risken H., Phys. Rev. A, 38 (1988) 2409. 

[10] Caldeira A. O. and Leggett A. J., Physica A, 121 (1983) 587. 

[11] Anastopoulos C. and Halliwell J. J., Phys. Rev. D, 51 (1995) 6870. 

[12] Zurek W. H. and Paz J. P., Phys. Rev. Lett, 72 (1994) 2508. 

[13] Ankerhold J., Europhys. Lett, 61 (2003) 301. 

[14] Isar A., Sandulescu A. and Scheid W., Int. J. Mod. Phys. B, 10 (1996) 2767. 

[15] Balescu R., Statistical dynamics (Imperial College Press, London) 1997. 



EUROPHYSICS LETTERS 



[16] Pathak A., J. Phys. A, 33 (2000) 5607. 
[17] Kandemir B. S., Phys. Lett. A, 245 (1998) 209. 

[18] Chen Y.-C. and Lebowitz J. L., Phys. Rev. B, 46 (1992) 10751; Phys. Rev. Lett, 69 (1992) 
3559. 



